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We study the effect of physicaf aging on the mechanicai properties of a model polymer glass 
using molecular dynamics simulations. The creep compliance is determined simultaneously with 
the structural relaxation under a constant uniaxial load below yield at constant temperature. The 
model successfully captures universal features found experimentally in polymer glasses, including 
signatures of mechanical rejuvenation. We analyze microscopic relaxation timescales and show 
that they exhibit the same aging characteristics as the macroscopic creep compliance. In addition, 
our model indicates that the entire distribution of relaxation times scales identically with age. 
Despite large changes in mobility, we observe comparatively little structural change except for a 
weak logarithmic increase in the degree of short-range order that may be correlated to an observed 
decrease in aging with increasing load. 
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I. INTRODUCTION 

Glassy materials are unable to reach equilibrium over 
typical experimental timescales [l|, H, Instead, the 
presence of disorder at temperatures below the glass tran- 
sition permits only a slow exploration of the configura- 
tional degrees of freedom. The resulting structural re- 
laxation, also known as physical aging is one of the 
hallmarks of glassy dynamics and leads to material prop- 
erties that depend on the wait time t w since the glass was 
formed. While thermodynamic variables such as energy 
and density typically evolve only logarithmically, the re- 
laxation times grow much more rapidly with wait time 

SIS- 

Aging is a process observed in many different glassy 
systems, including colloidal glasses @, microgel pastes 
0, and spin glasses 0, but is most frequently studied 
in polymers due to their good glass-forming ability and 
ubiquitous use in structural applications. Of particular 
interest is therefore to understand the effect of aging on 
their mechanical response during plastic deformation Q . 
In a classic series of experiments, Struik studied many 
different polymer glasses and determined that their stiff- 
ness universally increases with wait time. However, it has 
also been found that large mechanical stimuli can alter 
the intrinsic aging dynamics of a glass. Cases of both 
decreased aging (rejuvenation) Q and increased aging 
(overaging) [9J, [l0| have been observed, but the interpre- 
tation of these findings in terms of the structural evolu- 
tion remains controversial [ll], [l2j . 

The formulation of a comprehensive molecular model 
of the non-equilibrium dynamics of glasses has been im- 
peded by the fact that minimal structural change oc- 
curs during aging. Traditional interpretations of aging 
presume that structural relaxation is accompanied by a 



'Electronic address: mya@phas.ubc.ca 



decrease in free volume available to molecules and an as- 
sociated reduction in molecular mobility Q. While this 
idea is intuitive, it suffers from several limitations. First, 
the free volume has been notoriously difficult to define ex- 
perimentally. Also, this model does not seem compatible 
with the observed aging in glassy solids under constant 
volume conditions [X 31 ] . and cannot predict the aging be- 
havior under complex thermo-mechanical histories. Mod- 
ern energy landscape theories describe the aging process 
as a series of hops between progressively deeper traps in 
configuration space [lj,ll5|. These models have had some 
success in capturing experimental trends, but have yet to 
directly establish a connection between macroscopic ma- 
terial response and the underlying molecular level pro- 
cesses. Recent efforts to formulate a molecular theory of 
aging are promising but require knowledge of how local 
density fluctuations control the relaxation times in the 
glass jig ]. 

Molecular simulations using relatively simple models 
of glass forming solids, such as the binary Lennard- Jones 
glass [17] or the bead spring model [18] for polymers, 
have shown rich aging phenomenology. For instance, 
calculations of particle correlation functions have shown 
explicitly that the characteristic time scale for particle 
relaxations increases with wait time [191 ]. Recent work 
H3| has focused on the effect of aging on the mechan- 
ical properties; results showed that the shear yield stress 
(defined as the overshoot or maximum of the stress-strain 
curve) in deformation at constant strain rate generally in- 
creases logarithmically with t w . Based on a large number 
of simulations at different strain rates and temperatures, 
a phenomenological rate-state model was developed that 
describes the combined effect of rate and age on the shear 
yield stress for many temperatures below the glass tran- 
sition 

In contrast to the strain-controlled studies described 
above, experiments on aging typically impose a small, 
constant stress and measure the resulting creep as a func- 
tion of time and t w [H ■ In this study, we perform molecu- 
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lar dynamics simulations on a coarse grained, glass form- 
ing polymer model in order to investigate the relation- 
ship between macroscopic creep response and microscopic 
structure and dynamics. In Section IIII A[ we determine 
creep compliance curves for different temperatures and 
applied loads (in the sub-yield regime) and find that, as 
in experiments, curves for different ages can be super- 
imposed by rescaling time. The associated shift factors 
exhibit a power-law dependence on the wait time, and 
the effect of aging can be captured by an effective time 
as originally envisioned by Struik Q. In Section UlI Bl we 
compute microscopic mobilities and the full spectrum of 
relaxation times and show their relationship to the creep 
response. Additionally, we study several parameters that 
are sensitive to the degree of short-range order in Section 
IIII CI We detect very little evolution toward increased 
local order in our polymer model, indicating that short 
range order is not a sensitive measure of the mechanical 
relaxation times responsible for the creep compliance of 
glassy polymers. 



II. SIMULATIONS 

We perform molecular dynamics (MD) simulations 
with a well-known model polymer glass on the bead- 
spring level. The beads interact via a non-specific van 
der Waals interaction given by a 6-12 Lennard- Jones po- 
tential, and the covalent bonds are modeled with a stiff 
spring that prevents chain crossing [22]. This level of 
modeling does not include chemical specificity, but al- 
lows us to study longer aging times than a fully atomistic 
model and seems appropriate to examine a universal phe- 
nomenon found in all glassy polymers. All results will be 
given in units of the diameter a of the bead, the mass to, 
and the Lennard- Jones energy scale, u$. The characteris- 
tic timescale is therefore T£j = yjma 2 Juq, and the pres- 
sure and stress are in units of Uo/a 3 . The Lennard- Jones 
interaction is truncated at 1.5a and adjusted vertically 
for continuity. All polymers have a length of 100 beads, 
and unless otherwise noted, we analyze 870 polymers in a 
periodic simulation box. Results are obtained either with 
one large simulation containing the full number of poly- 
mers, or with several smaller simulations, each starting 
from a unique configuration, whose results are averaged. 
The large simulations and the averaged small simulations 
provide identical results. The small simulations are used 
to estimate uncertainties caused by the finite size of the 
simulation volume. 

To create the glass, we begin with a random dis- 
tribution of chains and relax in an ensemble at con- 
stant volume and at a melt temperature of 1.2uo/kB- 
Once the system is fully equilibrated, it is cooled over 
750txj to a temperature below the glass transition at 
T g w 0.35uo/fcs [HI]. The density of the melt is cho- 
sen such that after cooling the pressure is zero. We 
then switch to an NPT ensemble - the pressure and 
temperature are controlled via a Nose-Hoover thermo- 
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FIG. 1: Simulated creep compliance J(t,t w ) at a glassy tem- 
perature of T = 0.2tto/fcs for various wait times t w (indi- 
cated in the legend in units of txj). A uniaxial load of (a) 
a = OAuo/a 3 and (b) a = 0.5wo/a 3 is applied to the aged 
glasses. The strain during creep is monitored over time to 
give the creep compliance. 



stat/barostat - with zero pressure and age for various 
wait times (t w ) between 500 to 75,000tlj. The aged 
samples undergo a computer creep experiment where a 
uniaxial tensile stress (in the z-direction) is ramped up 
quickly over 75t£j, and then held constant at a value 
of cr, while the strain e = AL Z / L z is monitored. After 
an initial elastic deformation, the glass slowly elongates 
in the direction of applied stress due to structural relax- 
ations. In the two directions perpendicular to the applied 
stress, the pressure is maintained at zero. 



III. RESULTS 

A. Macroscopic Mechanical Deformation 

Historically, measurements of the creep compliance 
have been instrumental in probing the relaxation dynam- 
ics of glasses, and continue to be the preferred tool in 
investigating the aging of glassy polymers [ill HH, [24[ . In 
his seminal work on aging in polymer glasses, Struik Q 
performed an exhaustive set of creep experiments on dif- 
ferent materials, varying the temperature and the applied 
load. In this section, we perform a similar set of exper- 
iments with our model polymer glass. The macroscopic 
creep compliance is defined as 

J {t ,t w )= e -^-. (1) 
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Compliance curves J(t,t w ) for several temperatures and 
stresses were obtained as a function of wait time since 
the quench; representative data is shown in Figure Q] 
The curves for different wait times appear similar and 
agree qualitatively with experiment. An initially rapid 
rise in compliance crosses over into a slower, logarith- 
mic increase at long times. The crossover between the 
two regimes increases with increasing wait time. Struik 
showed that experimental creep compliance curves for 
different ages can be superimposed by rescaling the time 
variable by a shift factor, aj, 

J(t,t w ) = J(taj,t' w ). (2) 

This result is called time-aging time superposition 0, [j| . 
Simulated creep compliance curves from Fig. Q] can simi- 
larly be superimposed, and the resulting master curve is 
shown in Fig. [21 

Shift factors required for this data collapse are plot- 
ted versus the wait time in Fig. [3J All data fall along a 
straight line in the double- logarithmic plot, clearly indi- 
cating power law behavior: 

aj~t-». (3) 

This power law in the shift factor is characteristic of ag- 
ing, /z is called the aging exponent, and has been found 
experimentally to be close to unity for a wide variety of 
glasses in a temperature range near T g [J] . 

Figured] shows the effect of stress and temperature on 
the aging exponent, as determined from linear fits to the 
data in Fig. [31 At T = 0.2uo/kB, H is close to one for 
small stresses, but decreases strongly with stress. This 
apparent erasure of aging by large mechanical deforma- 
tions has been called "mechanical rejuvenation" [25[. Ex- 
periments have frequently found a stress dependence of 



FIG. 3: Plot of the shift factors found by superimposing the 
creep compliance curves, a.j (circles), the mean-squared dis- 
placement curves, clmsd (triangles), and the incoherent scat- 
tering function curves, ac (x) at different wait times (see 
text). The solid lines are linear fits to the data. 



the aging exponent [4| , although it is not always the case 
that the aging process slows down with applied stress; 
stress has been known to increase the rate of aging in 
some circumstances as well [j| The structural ori- 
gins of this effect are not well understood [Til, E2 • 

At T = 0.3uo/fcs, we find that the aging exponent is 
somewhat smaller than at T — 0.2uo/ks and varies much 
less with applied stress. This behavior is most likely 
due to the fact that the temperature is approaching T g . 
Indeed, experiments show that /j, rapidly drops to zero 
above T g . compliance is an order of magnitude larger at 
T = 0.3uo/ks than at T = 0.2uo//cs and the data does 
not fully superimpose in a master curve for long times 
where J > 0.2uo/a 3 . Shift factors were obtained from 
the small creep portion of the curves. 

The relatively simple relationship between shift factors 
and wait time permits construction of an expression that 
describes the entire master curve in Fig. [21 For creep 
times that are short compared to the wait time - such 
that minimal physical aging occurs over the timescale of 
the experiment - experimental creep compliance curves 
can be fit to a stretched exponential (typical of processes 
with a spectrum of relaxation times), 

J(t) = J cxp[(t/t ) m } (4) 

where io is the retardation factor, and the exponent, m, 
has been found to be close to 1/3 for most glasses 0|. 
A fit of this expression to our simulated creep compli- 
ance curves is shown in Fig. [2] (dotted line). This ex- 
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FIG. 4: The aging exponent, fi, determined from the slopes 
of log(aj) versus log(t w ) (from Fig. [3j plotted versus stress 
(open symbols). The solid symbols at zero stress refer to 
shift factors determined from clmsd (eq. [7)l and ac (eq. [6} 
data only. The dashed lines are guides to the eye. 



Struik's effective time, as it tries to subsume the modified 
aging dynamics during deformation in a single variable 
and in particular easily includes the effects of overaging 
or rejuvenation. 



B. Microscopic Dynamics 

The aging behavior of the simulated mechanical re- 
sponse functions agrees remarkably well with experiment. 
Additional microscopic information from simulations al- 
lows us to obtain more directly the relevant timescales 
of the system, and the relevant microscopic processes re- 
sponsible for aging. One parameter which has been use- 
ful in studying glassy dynamics is the "self part of the 
incoherent scattering factor [19], 

1 N 

C q (t,t w ) = —^2exp(iq- [?j{t w +t) - fj(t w )]) (6) 



pression is clearly only consistent with the data at times 
t < t w . At times much longer than the wait time, the 
creep compliance varies more slowly due to the stiffen- 
ing caused by aging during the course of the experiment. 
Struik suggested that eq. ([4]) could be extended to the 
long-time creep regime, where the experimental timescale 
may be longer than the wait time, by introducing an ef- 
fective time to account for the slowdown in the relaxation 
timescales: 

Upon replacing t with t c s, eq. (|4|) may be used to de- 
scribe the entire creep curve. Creep compliance curves 
from Fig. [2] can indeed be fit to this form (dashed lines) 
for a known wait time, t w , and aging exponent, fj,, as 
obtained from the master curve. We find m ~ 0.5 ±0.1 
for all stresses at T = 0.2ito/fc_B, and a relatively broad 
range of values for Jo and to are consistent with the data. 
For the simple thermo-mechanical history prescribed by 
the creep experiment, Struik's effective time formulation 
appears to work quite well. 

The present results parallel those of a recent simula- 
tion study of the shear yield stress in glassy solids [2~j| . 
In this work, the glassy solid was deformed at constant 
strain rate, and two different regimes of strong and weak 
rate dependence emerged depending on the time to reach 
the yield point relative to the wait time. In order to ra- 
tionalize these results, a rate-state model was developed 
that accounted for the internal evolution of the material 
with age through a single state variable 9(t). This formu- 
lation successfully collapses yield stress data for different 
ages and strain rates in a universal curve by adapting the 
evolution of the state variable during the strain interval. 
We note here that this state variable is closely related to 



where fj is the position of the j atom, and q is the 
wave-vector. C q curves as a function of age are shown 
in Fig. [5] and exhibit three distinct regions. Initially, 
C q decreases as particles make very small unconstrained 
excursions about their positions. There follows a long 
plateau, where the correlation function does not change 
considerably. In this regime, atoms are not free to dif- 
fuse, but are trapped in local cages formed by their near- 
est neighbours. For this reason, the time spent in the 
plateau regime is often associated with a "cage time" . 
The plateau region ends when particles finally escape 
from local cages (a- relaxation) , and larger atomic rear- 
rangements begin to take place. The cage time corre- 
sponds closely to the transition from short-time to long- 
time regime observed in the creep compliance. Structural 
rearrangements taking place in the a-relaxation regime 
are clearly associated with the continued aging observed 
in the creep compliance, as well as plastic deformations 
occurring in that region. 

The correlation functions for different ages are similar 
in form, but the time spent in the plateau region in- 
creases with age. Just as creep compliance curves can be 
shifted in time to form a master curve, we may overlap 
the long-time, cage-escape regions of C q by rescaling the 
time variable of the correlation data at different ages (see 
inset of Fig. O. The corresponding shift factors ac{t w ) 
are also shown in Fig. [3J where we see that the increase 
in cage time with age follows the same power law as the 
shift factors of the creep compliance. These results are 
qualitatively similar to the scaling of the relaxation times 
with age found in p^ ] with no load. 

The real space quantity corresponding to C q is the 
mean squared displacement, 

1 N 
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FIG. 5: Incoherent scattering factor (eq. [6} for different 
wait times measured under the same loading conditions as 
in Fig. DJa) for q = (0,0, 27r). The inset shows the master 
curve created by rescaling the time variable by ac- Symbols 
as in Fig. \V[h). 




FIG. 7: The displacement probability distribution versus time 
measured under the same loading conditions as in Fig. Ola) , 
with a wait time of 500tlj. The solid lines from left to right 
are obtained at times t of 75, 750, 7500, and 75000txj. The 
dashed lines show fits to the double Gaussian distribution (see 
text, eq. [8j) . 
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FIG. 6: Mean-squared displacement (eq. [TJl for different wait 
times measured under the same loading conditions as in 
Fig.QJa). The inset shows the master curve created by rescal- 
ing the time variable by clmsd- Symbols as in Fig. QJb). 



where Afj{t,t w ) — rj(t w + t) — fj(t w ). This function 
is shown in Fig. [6j Again we see three characteristic re- 
gions of unconstrained (ballistic) , caged, and cage-escape 
behavior. The departure from the cage plateau likewise 
increases with age, and a master curve can be constructed 
by shifting the curves with a factor cimsd (see inset of 
Fig. [6]). Shift factors clmsd are plotted in Fig. [3l along 
with shifts for creep compliance and incoherent scatter- 
ing function. As anticipated, the shifts versus wait time 
for (Ar 2 ) fully agree with those obtained from C q and J. 
This clearly demonstrates that for our model, the cage 
escape time is indeed the controlling factor in the aging 
dynamics of the mechanical response functions. 



Additional information about microscopic processes 
can be obtained by studying not only the mean of the dis- 
placements, but also the full spectrum of relaxation dy- 
namics as a function of time and wait time. To this end, 
we measure the probability distribution P(Ar(t, t w ) 2 ) of 
atomic displacements during time intervals, t, for glasses 
at various ages, t w . This quantity is complementary to 
the measurements of dynamical heterogeneities detailed 
in [26[ , where the spatial variations of the vibrational am- 
plitudes were measured at a snapshot in time to show the 
correlations of mobile particles in space. In our study, we 
omit the spatial information, but retain all of the dynam- 
ical information. 

Representative distribution functions are shown in 
Fig. [7] for a constant wait time of t w = 500t£j and var- 
ious time intervals t. The distributions were obtained 
from a smaller system of only 271 polymer chains due to 
memory constraints. The data does not reflect a simple 
Gaussian distribution, but clearly exhibits the presence 
of two distinct scales: there is a narrow distribution of 
caged particles and a wider distribution of particles that 
have escaped from their cages. This behavior can be de- 
scribed by the sum of two Gaussian peaks, 

/-Ar 2 \ /-Ar 2 \ 
P(Ar 2 ) = Ar lC xp^— —j+N 2 cxpl^— —J. (8) 

Deviations from purely Gaussian behavior are common in 
glassy systems and are a signature of dynamical hetero- 
geneities [HI, H3|- Experiments on colloidal glasses [HI 
show a similar separation of displacement distributions 
into fast and slow particles. 

A fit of the normalized distributions to eq. (J8j> (dashed 
lines in Fig. [7J requires adjustment of three parameters: 
the variance of caged and mobile particle distributions, 
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FIG. 8: The Gaussian fit parameters for the distribution of 
displacements (see text, eq. [8}, (a) Ni/N, (b) o\, and (c) erf 
measured under the same loading conditions as in Fig. [Ha) . 
The curves are for wait times increasing from left to right 
from 500t L j to 15000r LJ . 



<j\ and erf, as well as their relative contributions N\/N, 
where N = Ni + A^ . These parameters are sufficient to 
describe the full evolution of the displacement distribu- 
tion during aging. In Fig. [5J we show the fit parameters 
as a function of time and wait time. Again two distinct 
time scales are evident. At short times, most of the parti- 
cles are caged (Ni/N ~ 1), and the variance of the cage 
peak is also changing very little. There are few rear- 
rangements in this regime, however Fig. EJc) shows that 
a small fraction of particles are mobile at even the short- 
est times. At a time corresponding to the onset of cage 
escape, the number of particles in the cage peak begins 
to rapidly decay, and the variance of the cage peak in- 
creases. This indicates that the cage has become more 
malleable - small, persistent rearrangements occur lead- 
ing to eventual cage escape. In this regime, the variance 
of the mobile peak increases very little. Note that the 
typical length scale of rearrangements is less than a par- 
ticle diameter even in the cage escape regime, but the 
number of particles undergoing rearrangements changes 
by more than 50%. 

Similar to the compliance and mean-squared displace- 
ment curves, the data in Fig. [5ja) and (b) can also be 
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FIG. 9: The displacement probability distribution measured 
under the same loading conditions and wait times as in 
Fig. [lja) plotted at times corresponding to < r 2 (t,t w ) >= 
0.7, shown in the inset as a dashed line. The legend indicates 
the wait time. 



superimposed by shifting time. Shift factors for Ni/N 
and a\ coincide exactly with shifts for the mean; how- 
ever, data for erf (Fig. [3(c)) seems to be much less af- 
fected by the wait time. The aging dynamics appears to 
be entirely determined by the cage dynamics, and not by 
larger rearrangements within the glass. 

Since the fit parameters exhibit the same scaling with 
wait time as the mean, one might expect that the en- 
tire distribution of displacements under finite load scales 
with the evolution of the mean. In Fig. we plot dis- 
placement distributions for several wait times at time in- 
tervals chosen such that the mean squared displacements 
are identical (see inset). Indeed, we find that all curves 
overlap, indicating that the entire relaxation spectrum 
ages in the same way. A similar observation was recently 
made in simulations of a model for a metallic glass aging 
at zero stress [2!| , although in this study the tails of the 
distribution were better described by stretched exponen- 
tials. 

In order to study the effect of load on the relaxation 
dynamics, we compare in Figure [TU] the fit parameters for 
a sample undergoing creep (replotted from Fig. [5]) and a 
reference sample without load. It is clear that the dy- 
namics are strongly affected by the applied stress only in 
the region characterized by a-relaxations. For the stress 
applied here, the onset of cage-escape does not appear to 
be greatly modified by the stress, however the decay in 
Ni/N and the widening of the cage peak are accelerated. 
The stress does not modify the variance of the mobile 
peak, confirming again the importance of local rearrange- 
ments as compared to large-scale motion in the dynamics 
of the system. The accelerated structural rearrangements 
caused by the stress result in creep on the macroscopic 
scale, but may also be responsible for the modification of 
the aging dynamics with stress as observed in Fig. 2J 
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FIG. 10: The Gaussian fit parameters to the displacement 
distributions (see text, eq. [8| (a) Ni/N, (b) af, and (c) a\ of 
a sample aged at T — 0.2uo/fcs for t w = 500tlj, and then 
either undergoing a creep experiment at a = 0.4uo / a 3 (black) , 
or simply aging further at zero stress (red). 



C. Structural evolution 

The connection between the dynamics and the struc- 
ture of a glass during aging remains uncertain, mostly 
because no structural parameter has been found that 
strongly depends on wait time. Recent simulation stud- 
ies of metallic glasses have shown the existence of sev- 
eral short range order parameters that can distinguish 
between glass y st ates created through different quench- 
ing protocols [30l . [3ll l32l | . A strong correlation has been 
found between "ordered" regions of the glass and strain 
localization. Many metallic glasses are known to form 
quasi-crystalline structures that optimize local packing. 
It remains to be seen whether the short-range order 
evolves in the context of aging and in other glass for- 
mers such as polymers and colloids. A recent experimen- 
tal study of aging in colloidal glasses found no change in 
the distribution function of a tetrahedral order parame- 
ter [33[ . Below, we investigate several measures of local 
order in our model as they evolve with age and under 
load. 

Since Lennard- Jones liquids are known to condense 
into a crystal with fee symmetry at low temperatures, 
it is reasonable to look for the degree of local fee order in 



our polymer model. The level of fee order can be quan- 
tified via the bond orientational parameter (34|, 

Q« = (if E M 2 ) • P) 

This parameter has been successfully used to character- 
ize the degree of order in systems of hard sphere glasses. 
Qe is determined for each atom by projecting the bond 
angles of the nearest neighbours onto the spherical har- 
monics, Ye m (9, (f>). The overbar denotes an average over 
all bonds. Nearest neighbours are defined as all atoms 
within a cutoff radius, r c , of the central atom. For all of 
the order parameters discussed here, the cutoff radius is 
defined by the first minimum in the pair correlation func- 
tion, in this case 1.45a. Qg is approximately 0.575 for a 
perfect fee crystal; for jammed structures, it can exhibit 
a large range of values less than about 0.37 [34|]. The full 
distribution of Qq for our model glass is shown for several 
ages as well as an initial melt state in Fig. [TTT a). We see 
that there is very little difference even between melt and 
glassy states in our model, and no discernible difference 
at all with increasing age. 

Locally, close-packing is achieved by tetrahedral order- 
ing and not fee ordering, however, tetrahedral orderings 
cannot span the system. The glass formation process has 
been described in terms of frustration between optimal 
local and global close-packing structures. To investigate 
the type of local ordering in this model, we investigate 
a 3-body angular correlation function, P(ff). 9 is de- 
fined as the internal angle created by a central atom and 
individual pairs of nearest-neighbours, and P(9) is the 
probability of occurrence of 9. Results for this corre- 
lation are shown in Fig. Illf b). Two peaks at approxi- 
mately 60° and 110° indicate tetrahedral ordering. The 
peaks sharpen under quenching from the melt, but the 
distribution does not evolve significantly during aging. 
In contrast, simulations of metallic glass formers showed 
a stronger sensitivity of this parameter to the quench 
protocol [3l|, but most of those changes may be due to 
rearrangements in the supercooled liquid state and not 
in the glassy state. 

Another parameter that has been successful in classi- 
fying glasses is the triangulated surface order parameter 

s = E( 6 -<?K (10) 

Q 

which measures the degree of quasi-crystalline order. The 
surface coordination number, q, is defined for each atom 
of the coordination shell as the number of neighbouring 
atoms also residing in the coordination shell; v q is the 
number of atoms in the coordination shell with surface 
coordination q. Ordered systems have been identified 
with S equal to 12 (icosahedron) , 14, 15 and 16. Figure 
\TT\ c) shows the probability distribution for P(S) for the 
melt and for glassy states with short and long wait times. 
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the nature of the structural changes during rejuvenation. 




FIG. 12: Precent change in the triangulated surface order 
parameter with wait time as compared to the just-quenched 
sample. Circles are for samples aged at zero pressure for the 
time t w . Triangles are for the same samples immediately after 
ramping up to the creep stress. For T = 0.2uo/fcs this stress 
is OAuo/a 3 , and for T = 0.3uo/fcs the stress is O.lwo/a 3 . 
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FIG. 11: Short-range order parameters: (a) the bond- 
orientational parameter, (b) the three-body angular corre- 
lations, (c) the surface triangulated order (see text for dis- 
cussion), x's show the melt state, circles show the sample 
aged for t w = 500tlj, and triangles show the sample aged for 
500,000t lj . 



The peak of the distribution moves toward lower S (more 
ordered) upon cooling, and continues to evolve slowly in 
the glass. The mean of S relative to the as-quenched 
state, (S), is shown in Fig. as a function of wait time 
at two temperatures. We see that (S) is a logarithmi- 
cally decreasing function of wait time. Even though this 
is not a strong dependence, this order parameter is sig- 
nificantly more sensitive to age than the others that have 
been investigated. 

Figure H21 also shows the order parameter (S) after the 
ramped- up stress has been applied to the aged samples. 
We can see that at T — 0.2ito/fes, some of the order 
that developed during age is erased, while no appreciable 
change occurs at the higher temperature T = 0.3wo/fcs. 
These results correlate well with the behavior of the ag- 
ing exponent found in Fig. 31 where mechanical rejuve- 
nation was found at lower temperatures but was much 
less pronounced at higher T. The load is applied very 
quickly, and most of the deformation in this regime is 
affine, however, the strain during this time was similar 
for both temperatures, therefore the effect is not simply 
due to a change in density. More work is needed to clarify 



IV. CONCLUSIONS 

We investigate the effects of aging on both macroscopic 
creep response and underlying microscopic structure and 
dynamics through simulations on a simple model polymer 
glass. The model qualitatively reproduces key experi- 
mental trends in the mechanical behavior of glasses under 
sustained stress. We observe a power-law dependence of 
the relaxation time on the wait time with an aging expo- 
nent of approximately unity, and a decrease in the aging 
exponent with increasing load that indicates the presence 
of mechanical rejuvenation. The model creep compliance 
curves can be fit in the short and long-time regimes using 
Struik's effective time formulation. Additionally, inves- 
tigation of the microscopic dynamics through two-time 
correlation functions has shown that, for our model glass, 
the aging dynamics of the creep compliance exactly corre- 
sponds to the increase in the cage escape or a-relaxation 
time. 

A detailed study of the entire distribution of parti- 
cle displacements yields an interesting picture of the mi- 
croscopic dynamics during aging. The distribution can 
be described by the sum of two Gaussians, reflecting 
the presence of caged and mobile particles. The frac- 
tion of mobile particles and the amplitude of rearrange- 
ments in the cage strongly increase at the cage escape 
time. However, in analogy with results in colloidal glasses 
[HI , structural rearrangements occur even for times well 
within the "caged" regime, and fairly independent of wait 
time and stress. For our model glass, we find that the 
entire distribution of displacements scales with age in 
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the same way as the mean. At times when the long- 
time portion of the mean squared displacement overlaps, 
the distribution of displacements at different wait times 
completely superimpose, confirming that all of the me- 
chanical relaxation times scale in the same way with age. 

To characterize the evolution of the structure during 
aging, we investigate several measures of short-range or- 
der in our model glass. We find that the short-range order 
does not evolve strongly during aging. The triangulated 
surface order [32| . however, shows a weak logarithmic de- 
pendence on age. Results also show a change in structure 
when a load is rapidly applied, and this seems to be cor- 
related with an observed decrease in the aging exponent 
under stress. 

This study has characterized the dynamics of a model 
glass prepared by a rapid quench below T g , followed by 
aging at constant T and subsequent application of a con- 
stant load. For such simple thermo-mechanical histo- 
ries, existing phenomenological models work well, how- 
ever, the dynamics of glasses are in general much more 



complex. For instance, large stresses in the non-linear 
regime modify the aging dynamics and cause nontrivial 
effects such as mechanical rejuvenation and over-aging 
[ToL [ll[ . Also, experiments have shown that the time-age 
time superposition no longer holds when polymer glasses 
undergo more complex thermal histories such as a tem- 
perature jump [23TJ . The success of our study in analyzing 
simple aging situations indicates that the present simula- 
tion methodology will be able to shed more light on these 
topics in the near future. 
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